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INTRODUCTION 

A numerical method was developed by Dey (1977) 
for solving nonlinear systems, some applica- 
tions of which were later made to stiff sys- 
tems (Dey, 1982). Convergence analysis was 
done using nonlinear D-mappings (Dey, 1981). 


G k : D k+i * D k D k+i • 

G (u ,u ) - G k (u*,u*) 

” A k (u k+1 - u*) + B k (u k _ u*) (3) 


It is extremely difficult to represent this 
analysis computationally. Local linearization 
for such an analysis, which rendered computa- 
tional modeling of D-mappings feasible, was 
suggested by Lomax (1983). In this article 
we discuss linearized modeling of D-mappings 
and some applications of the method. 


D-MATRICES AND D-MAPPINGS 

If a sequence of square matrices of the same 
order satisfy the following condition, 

lim Vk-i • • ' A i = 0 (1 > 

k-*-« 


each A k is called a D-matrix. A D-matrix 
is not necessarily a convergent matrix, and 
conversely. 


and V k » K, and if (I - A k ) _1 B k is a 
D-matrix, G k is called a D-mapping (Dey, 
1981). 

If we now consider a nonstationary iterative 
scheme of the form 


k+i 


„ . k+i k. 
G k (u ,u ) 


(4) 


and if G k . D k+1 * D k ■* D k+1 is a D-mapping, 
then 

lim u k “ u* (5) 

k-*=> 

where u* = G k (u*,u*) V k (Dey, 1981). 


LOCALLY LINEARIZED D-MAPPING 


Theorem 1 . A sufficient condition that A k 
is a D-matrix is that 

ll\ !l q S a < 1 (2) 

V k > K and chat q is the same V k. 


Let us linearize (4) on D k 

first-order approximation of 

near (u k ,u k ). Then, 

k+2 _ , k k. . 

u = G k^ u ,u + G k^ u 


x D k , using 

r , k+i 
G k (u ,u 

k+2 k. 

- u ) 


( 6 ) 


Theorem 1 is easily proved. Let 


k , k k 
“ = (U 2 U 2 


U J> T £ “k 


k = 1, 2, . ., (uj = value of uj at some 

kth iteration). 


Let us consider a chained linear spaces 

D k £ °k -2 C . . . c D C R n . R n = n-dimensional 

real space. Let u k , u* € D k V k and 


where G k is the Frdchet derivative of G k 
on D k x £> k . Equation (6) may be expressed as 

k+2 . k , , ,,, 

u = Aj^u + b k (7) 

where A k = -(I - G^)'^, 

b k = (I - G k ) _1 G k (u k ,u k ) . We have assumed 
that (I - G k ) is invertible. Now we may 
prove a second theorem. 
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Theorem 2 . If (1) |Aj - A*| < E, where E 
is a matrix consisting of elements that are 
positive and arbitrarily small and 
(ii) jbj - b*| < e, e is a vector consisting 
of elements that are positive and arbitrarily 
small, then (5) is true (convergence) if A^ 
is a D-matrix (Dey, 1983a). 

Theorem 3 . If Gfe is a D-matrix, so is A^ 
(Dey, 1983a). 

This principle may now be applied 
computationally. 

PERTURBED FUNCTIONAL ITERATION 

Let a nonlinear system be expressed as 


Following Theorems 2 and 3 we may prove that 
if is a D-matrix, then (12) is true 

(linearized sense). Recent results (Dey, 
1983b) using local linearization indicates 
that if 

' max|G | s B/J (13) 

j ,m J 

where Gj m = 3Gj /Sujj and 0 < 8 < 1, G in 
(9) is a D-mapping. In order that (13) may 
be correct, certain input parameters for the 
system (e.g., mesh size and time-step) have 
to be chosen m special ways. If this cannot 
be found a convex-type operation may be 
defined as follows: 


u = G 0 (u) 
u € D, G 


( 8 ) 


A Gauss-Seidel-type iteration for the solu- 
tion may be expressed as 


k+i 


k+i k. 
G(u ,u ) 


(9) 


G: D x d -*■ D, u K £ D V k 


A perturbed iterative scheme (Dey, 1977) may 
be expressed as (in the element form) 


k+i k+i , k+l k, 

i • o). + G. (u , u ) 

J J 3 


( 10 ) 


where 


k+i 

"J 


[G j (Gj +1,k ) 


_k+i ,k. 

G j 1 


x (1 


3 G k+i,k r i 
J J 


(Ha) 


_k+i,k _ , k+i k+i k 

Gj “Vu, • * • Vi* tt J 


kv 

• u j> 

(lib) 


G j(Gj +1 ' k ) = G j (u ^ +1 


k+i 


_k+i ,k k 


j 


■Uj+, • • • U J> 


(11c) 


v G r* k) - t3G j /3u 3 ] k+i k+c 

u ! ,u 2 • • • 

,,k+i r k+i,k 

V l ’l 

k k 

u j+i ' • - U J 

(lid) 

The oij term is a perturbation parameter 
which accelerates the rate of convergence of 
(9) and stabilizes the numerical algorithm. 

It has been proved (Dey, 1981) that if G is 
a D-mapping on Dj c+1 x D^, a necessary and 
sufficient condition for convergence is 


lim a) 
k~ J 


0 V j 


( 12 ) 


G^ (u,u) = (1 - ct^ )u^ 


+ a G 
3 i 


(u,u) 


(14) 


Assuming Gjj(u,u) j* 1, it has been found 
that G is a D-mapping (locally linearized) 
for the following: 


1. oj = (-l)P(l - Gj ) 1 if Gj m = 0, 

m i j and p = 0 if G.. < 1, p ■ 1 if 

Gj j > !• 


2. a = min 
J l<m<j 
ra^j 


/ 8/J , ,,p (1 + 8/J) i 

VV • 1 - h J 


if Gj m i 0, mi< j and p = 0 if G.. < 1 

p = 1 if Gj j > 1. 


where B is such that 

1 1 - G. 


I > 6 > J 1 + 


-U_ 


jm I 


(15a) 


(15b) 


3. Oj = 1 if (13) is satisfied. 

The algorithm of perturbed functional itera- 
tion including a linearized convergence 
analysis may be briefly expressed as follows. 
At each iteration level, compute Gj m , 
m = 1, 2, . . .,J. If (13) is satisfied, set 
a j = 1 ; otherwise, compute Oj using (15). 

If Oj # 1, replace Gj by Gj , as given by 
(14). Compute ti) j using (11a) — (lid) and 
compute u, at the new iteration level by 
(10). If (12) is satisfied at some iteration 
level, convergence is found; l f G 1J " 1 , the 
method fails. 


In general, for a J*J system the method 
requires (i) J‘ + J functionals to be com- 
puted for convergence analysis, (ii) partial 
linearization along the diagonal, and 
(lii) no Jacobians. 

It has been proved analytically (Dey, 1977) 
that m the vicinity of the root, the method 
should display a superlinear rate of 
convergence . 
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A DEGENERATE IMPLICIT CODE 
Let a nonlinear model be represented by 

du/dt = f(u), u = (u 2 u 2 . . . u/ 

(16) 

u(0) = u Q (initial condition). Approximating 
(16) by a two-point backward-difference 
scheme, we get: 

n+i n . . . . n+i. , , 

Uj = Uj + Atfj(u ) , At = time-step 

(17) 

This nonlinear system may now be solved by 
the above method which forms a degenerate 
implicit code (since the one-step, matrix- 
inversion principle is not used for solution). 
If a convex-type operation of the form (14) 
is used, (16) becomes 

du/dt = (I - a)du/dt + af(u) (18) 

where a “ diaglaj.Oj, . . . oj) and 1 = the 
identity matrix. If (16) is a stiff system, 
(condition number of f'(u) >> 1), it 
generally requires At to be very small if a 
functional iteration of the form (9) is 
applied for solution. In (18), a scales the 
elements in the Jacobian matrix f'(u), and 
using (15) for a. 's mean (for a given At), 
D-mapping is found so that perturbed func- 
tional iteration converges. 


APPLICATIONS 
Application 1 . 

(Bui, 1979) Uj = -10,004 Ul + 10,000 u 2 ; 
u 2 = Uj - u 2 - uj; u j (0) = u 2 (0) = 1. An 
approximate solution is 
u 2 (t) = (10,004 exp(-3t) / [ 10,008 
- 4 exp(-3t) ] } 1 / 3 

Uj(t) = (10, 000/10, 004)uJ . As t -*■*>, u L (t) -*■ 0, 
u 2 (t) -*■ 0. Linearizing this system near 
u 1 (0),u 2 (0), the condition number is 10** . 

Using linearized convergence analysis for the 
degenerate implicit code we got At z 10” 5 , 
a sufficient condition for convergence if 
a-j = 1. Introducing (18) and computing cij ' s 
in a subroutine using (15) we used At = 10 -3 
to 1O 0 ; correct results were found. No pro- 
gram interruption was cased. Details may be 
found in Dey (1983b). 

Application 2. Irradiation of Neutral Water 

The model developed by Chatterjee and Magee 
and che analysis of its numerical solution 
are given m Chatterjee and others (1983). 

The equations and the rate constants are 
given in Table 1. For our present analysis 
we linearized the system and computed aj's. 
Stiffness was measured by Strate (1983) at 
c = 0, 0.1, 1, and 10. Condition numbers are, 
respectively, 10 13 , 10 12 , 10 12 , ^^(approxi- 
mately). This may be seen to be true in 
Fig. 1. This pattern of solution was ana- 


lyzed by Chatterjee and Magee and was found 
to be valid. Here, difference equations weie 
formed by approximating the derivatives by 
using the two-point trapezoidal rule. 
D-mappings were introduced, and time-accurate 
solutions were computed with At = 10 -8 , 10 -6 . 


CONCLUSION 

Numerical solutions of stiff systems are 
generally obtained by using multistep 
implicit codes (Miranker, 1981) which require 
inversion of matrices obtained by computing 
Jacobians. This has been avoided in the 
technique explained here. However, the code 
is dependent on the Jacobians for its conver- 
gence analysis Such a linearized analysis 
seems to be quite effective, and, in con- 
trast with its nonlinear counterpart, the 
complete analysis can be done computation- 
ally. More applications are under 
consideration. 
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Table 1 Differential Equations for Transient Species and 
Radiation Products in Irradiated Neutral Water 


^p. = ^(H ) 2 - k 2 (e' Aq )(H) - k s (H) (OH) + k 7 (H 3 0 + ) (e“ ) - k 9 (H) (H 2 0 2 ) 

+ k 12 (OH)(H 2 ) - kj 3 (HOj ) (H) - k 17 (H)(0 2 ) 

(ft ( e_ Aq) “ 'M e ' A q )(H) - 2k 3 (e ‘ Aq ) 2 - k,(e' Aq )(OH) - k 7 (H 3 0) (e' Aq ) 

- kio<e' Aq >(H 2 0 2 ) - k l4 (e‘ Aq )(0 2 ) 

(OH) = -k 4 (e~ Aq ) (OH) - k 5 (H) (OH) - 2k 8 (OH ) 2 + k,(H)(H 2 0 2 ) + k 2 „ (e" A ) (H 2 0 2 ) 

- k 2 j (OH) (HjjOj ) - k 12 (0H)(H 2 ) - k l 5 (H0 2 )(0H) 

(H>) = -k 7 (H 3 0)(e" Aq ) - k 3 (H 3 0 + )(OH‘) - k la (H 3 0)(0 2 ") + k 19 (H 2 0) V k 2 f) (H0 2 ) 

^ (H z O) = k 5 (H) (OH) + k 8 (H 3 0)(0H') + k 9 (H)(H 2 0,) + k 2 2 (OH) (H 2 0 2 ) + k J 2 (OH)(H 2 ) 
+ k 15 (H0 2 )(0H) + k 18 (H 3 0)(0 2 ') - k 19 (H 2 0) 

(H 2 ) « kj (H ) 2 + k 2 (e' Aq )(H) + k 3 (e ' Aq ) 2 - k 12 (OH)(H 2 ) 

^(H;0 2 ) - kg (OH ) 2 - k,(H)(H 2 0 2 ) - k 10 (e' Aq )H 2 O 2 ) - k 3 3 (OH) (H,0 2 ) + k 13 (H0 2 )(H) 
+ k 16 (H0 2 ) 2 


+ 0.55 I 
+ 2.65 I 

+ 2.70 I 
+ 2 65 I 

- 4.10 I 
+ 0.45 1 

+ 0.70 I 


(OH') = k 2 (e~ Aq )(H) + 2k 3 (e ~ Aq ) 2 + k„(e A )(OH) - k 8 (H 3 0)(0H') 
+ k 10 (e' Aq )(H 2 O 2 ) + k 19 (H 2 0) 


h < H0 *> 

« k ll (0H)(H 2 0 2 ) - 

k 13 (H0 2 )(H) 

- k 15 (H0 2 )(0H) - 2k 18 (H0 2 ) 2 

+ k 17 

(H)(0 2 ) 

+ k l 8 (H 3 0) (0 2 ~) 

“ ^ 20 (^ 2 ) 




Tt 

* ' k i- (e ' Aq )( 0 2 ) 

+ k l 5 (H0 2 )(OH) + k 16 (H0 2 ) 2 - k 17 (H)(0 2 ) 



h ( ° 2_) 

= k l4 (e Aq ) ( 0 2 ) 

k l e (H 3 O) (o 2 

") + k ; 0 (HO,) 



where 

k 3 = 10 10 

k 2 13 

2.5xl0 10 k 3 = 6 X 10 9 

k . 4 

= 3 X 10 


k s = 2.4xl0 10 k 8 = 

4 X 10 9 k 7 = 2 . 3 X 10 1 0 

kg 

= 3 X 10 


k 9 = 10 3 

k io = 

1 . 2 X I0 l 0 k 13 = 5 X 10 7 

k l 2 

= 6 X 10 


k i3 = 10l ° 

k 14 = 

1 . 9 X 1 0 1 0 k 15 = 10 10 

k l 6 

= 2 x 10 ' 


k 1 7 = 10 10 

k i 0 

3 X 10 10 k , 9 = 5 5 X 10 -5 

k 2 o 

= 10 6 
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TIME, sec 

Fig. 1 Concentrations of species (Ex 2) vs. time for I = 6.667 < 10 -7 in the 
logarithmic scale up to t = 30 sec. (Here steady state is reached for 
all the species.) 
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